NWS Forecast Zones

Summary

This notebook imports the National Weather Service Forecast Zones, as potential spatial units for defining groups of CIMIS stations for data cleaning.

The NWS issues forecasts and some watches and warnings for public zones which usually are the same as counties but in many cases are subsets of counties. Counties are subset into zones to allow for more accurate forecasts because of the differences in weather within a county due to such things as elevation or proximity to large bodies of water.

Setup

Load a few libraries:

library(sf)
library(dplyr)
library(tigris)
options(tigris_use_cache = TRUE)
library(cimir)
library(tmap)

Import the Forecast Zones

The Forecast Zones are available as a zipped up Shapefile. We can use gdal’s virtual file system prefixes to download and unzip it.

publiczones_zip_url <- "https://www.weather.gov/source/gis/Shapefiles/WSOM/z_05mr24.zip"
publiczones_vfs_url <- paste0("/vsizip", "/vsicurl/", publiczones_zip_url)
publiczones_vfs_url

[1] "/vsizip/vsicurl/https://www.weather.gov/source/gis/Shapefiles/WSOM/z_05mr24.zip"

publiczones_sf <- st_read(publiczones_vfs_url, 
                          query = "SELECT * FROM z_05mr24 WHERE STATE = 'CA'") 

Reading query `SELECT * FROM z_05mr24 WHERE STATE = 'CA''
from data source `/vsizip/vsicurl/https://www.weather.gov/source/gis/Shapefiles/WSOM/z_05mr24.zip' 
  using driver `ESRI Shapefile'
Simple feature collection with 173 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.4096 ymin: 32.53566 xmax: -114.1312 ymax: 42.00951
Geodetic CRS:  NAD83

head(publiczones_sf)

Import County Boundaries

ca_counties_sf <- tigris::counties(state = "CA", cb = TRUE, resolution = "5m") |> 
  select(GEOID, NAME)

Retrieving data for the year 2022

ca_counties_sf

Import Active CIMIS Stations

cimis_stn_rds_fn <- "cimis_stn.Rds"

if (file.exists(cimis_stn_rds_fn)) {
  cimis_stn_sf <- readRDS(cimis_stn_rds_fn)

} else {
  
  cimir::set_key("abcd")
  cimis_stn_sf <- cimir::cimis_station() |> 
    filter(IsActive == "True") |> 
    select(stid = StationNbr, name = Name, HmsLatitude, HmsLongitude) |> 
    distinct() |> 
    mutate(stid = as.numeric(stid),
           lon = as.numeric(gsub("^.*/ ", "", HmsLongitude)),
           lat = as.numeric(gsub("^.*/ ", "", HmsLatitude))) |> 
    select(stid, name, lon, lat) |> 
    sf::st_as_sf(coords = c("lon", "lat"), crs = 4269) 

  saveRDS(cimis_stn_sf, file = cimis_stn_rds_fn)
  
}

Overlay Forecast Zones, CIMIS Stations, and County Boundaries

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(publiczones_sf |> transmute(name_zone = paste0("NWS Zone: ", NAME, ", ", ZONE))) +
  tm_polygons(col = "MAP_COLORS", alpha = 0.5) +
tm_shape(ca_counties_sf |> select(NAME, GEOID)) +
  tm_borders(col = "#555", lwd = 2) +
tm_shape(cimis_stn_sf |> transmute(stid_name = paste0("CIMIS #", stid, ": ", name))) +
  tm_dots(col = "blue", clustering = FALSE) 

Conclusions

Although I don’t know any details of how the NWS constructed Forecast Zones, from the description they sound like they could be candidates for grouping CIMIS stations together.

The next steps might be to overlay the CIMIS stations on the Forecast Zones and see how they line up.